arXiv:1507.05078v2 [cond-mat.quant-gas] 4 Nov 2015 


A Non-Equilibrium Kinetic Theory for Trapped Binary Condensates 


M. J. Edmonds, K. L. Lee, and N. P. Proukakis 
Joint Quantum Centre (JQC) Durham-Newcastle, School of Mathematics and Statistics, 
Newcastle University, Newcastle upon Tyne NE1 7 RU, England, UK 
(Dated: November 5, 2015) 

We derive a non-equilibrium finite-temperature kinetic theory for a binary mixture of two interacting 
atomic Bose-Einstein condensates and use it to explore the degree of hydrodynamicity attainable in 
realistic experimental geometries. Based on the standard separation of timescale argument of kinetic 
theory, the dynamics of the condensates of the multi-component system are shown to be described by 
dissipative Gross-Pitaevskii equations, self-consistently coupled to corresponding Quantum Boltz¬ 
mann equations for the non-condensate atoms: on top of the usual mean field contributions, our 
scheme identifies a total of 8 distinct collisional processes, whose dynamical interplay is expected 
to be responsible for the system’s equilibration. In order to provide their first characterization, we 
perform a detailed numerical analysis of the role of trap frequency and geometry on collisional rates 
for experimentally accessible mixtures of 8l Rb- 41 K and s 'Rb- 85 Rb, discussing the extent to which 
the system may approach the hydrodynamic regime with regard to some of those processes as a 
guide for future experimental investigations of ultracold Bose gas mixtures. 
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I. INTRODUCTION 

The possibility of unprecedented control over ex¬ 
perimental parameters in ultracold atom experiments, 
such as the statistics, interactions and dimensional¬ 
ity of trapped gases m, offers an opportunity to 
elucidate novel many-body quantum effects. At the 
level of a single-component Bose gas, the study of the 
Bose-Einstein condensate has already bifurcated into a 
plethora of directions. Opportunities now exist to in¬ 
vestigate a broad spectrum of problems, including mim¬ 
icking the behaviour of atoms in solids using sophisti¬ 
cated optical manipulations m, as well as applica¬ 
tions to quantum information and computation [SMI- 
Experimental advances have also led to the controlled 
generation of multi-component |11H23] and spinor [24J- 
[50] condensates, with the dynamical interplay between 
different components leading to even richer physics, in¬ 
cluding, for example, phase separation [SlTfoSI and spin- 
domain formation [25] 29.. Recently, condensates have 
also been used to simulate gauge theories, which has at¬ 
tracted intense experimental and theoretical focus due 
to the strong analogies with condensed matter systems 
[34 . 351 . The behavior of Bose gas mixtures is also re¬ 
lated to the study of doubly-superfluid Bose-Fermi mix¬ 
tures in the BEC regime 3f>‘, where the Fermi gas forms 
a molecular condensate. 

For the single-component condensates, an understand¬ 
ing of the dynamics of Bose-condensed systems often re¬ 
lies on the Gross-Pitaevskii equation, which naturally en¬ 
compasses the wave-like behaviour of the weakly inter¬ 
acting gas, valid deep within the ultracold regime. How¬ 
ever, to gain insight into the dynamics of the gas over a 
broader range of temperatures, one must explicitly con¬ 
sider the behaviour of the normal component of the sys¬ 
tem, which leads to a rich non-equilibrium behaviour. 
Numerous approaches have been devised to describe the 


condensate dynamics in the presence of a thermal cloud 
(see e.g the reviews EIH12), each with its own merits 
and drawbacks. Classical field methods [4lT - [46] cumula¬ 
tively describe the highly-occupied low-lying ’’classical” 
inodes of the gas, relying on the ergodic relaxation of a 
non-equilibrium initial state (to a Rayleigh-Jeans distri¬ 
bution) ; appropriately-sampled quantum noise could also 
be added in the initial states to mimic quantum fluctu¬ 
ations (an approach referred to as ’’truncated Wigner”) 
[miii]. Explicitly adding a stochastic coupling to a heat 
bath, representing the set of high-lying modes largely un¬ 
affected by the condensate, one can also introduce fluctu¬ 
ating dynamics into the system [45H55] ; this is expected 
to be mostly relevant for studying equilibrium fluctua¬ 
tions [54U59] and quenched dynamics j52[ ffiOHfiJ] . While 
such approaches are suited for describing the critical re¬ 
gion, they only describe dynamics up to a (fixed en¬ 
ergy/momentum) cutoff 141] . and cannot therefore ac¬ 
count for any perturbations of the high-lying, thermal, 
inodes. 

Contrary to such approaches, the dynamics of ther¬ 
mal modes can be accurately handled by an alternative 
perturbative method, following the usual route of kinetic 
theory, which describes the coupled condensate and ther¬ 
mal cloud dynamics, based on a separation of timescales 
argument |641l73j ; while this method relies on symmetry- 
breaking m , and thus fails to account for the criti¬ 
cal fluctuation region, it is particularly suited to study¬ 
ing damping of collective modes and macroscopic excita¬ 
tions, which it has done very successfully [75M81] . De¬ 
spite its inherent limitation in requiring the assumption 
of a non-zero condensate mean field (which can however 
be negligibly small), this method (referred to by many 
as the ” Zaremba-Nikuni-Griffin”, or ZNG method m) 
has nonetheless been found to perform very well even 
on the issue of condensate number growth following a 
sudden truncation in the thermal distribution [82] or on 
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surface evaporative cooling [83], complementing studies 
based on other approaches which do not themselves re¬ 
quire a symmetry-broken condensate mean field potential 
when initiating the numerical simulations j35J [BIJ I51H87] . 

A somewhat similar kinetic approach, which is explic¬ 
itly number-conserving, and does not invoke symmetry¬ 
breaking [88] has also been successfully implemented for 
describing system dynamics HHHg. 

In the context of multi-component condensates, which 
have been extensively studied with coupled Gross- 
Pitaevskii equations (GPEs) [5IH331 . or their dis¬ 

sipative generalisations [1001 - 1102] . their finite temper¬ 
ature dynamics remains a partly open problem. Ap¬ 
proaches considered to date include classical field [103], 
truncated Wigner [104H106] . coupled stochastic pro¬ 
jected Gross-Pitaevskii equations [1071 - 11101 . or number- 
conserving approaches m- However, the detailed dy¬ 
namics far from the critical region are expected to be 
better described by a model that fully accounts for all 
condensate and thermal cloud dynamics. This is par¬ 
ticularly important since, parallel to the internal relax¬ 
ation within each system, the two dynamical thermal 
clouds will also need to equilibrate together, thus creat¬ 
ing a rather involved competition of collisional processes, 
with distinct timescales. While the promising number- 
conserving method of Ref. m has not yet been ad¬ 
vanced to the self-consistent dynamical level, all other 
methods (classical field, truncated Wigner and stochas¬ 
tic GPEs) feature a cutoff, and thus ignore the coupling 
of the high-lying thermal modes within and across the 
two systems; although such an approximation may be 
adequate for certain non-equilibrium features (e.g. de¬ 
fect formation following a quench [60] , persistent current 
decay [112] 1. it is nonetheless known to fail, at least for¬ 
mally, in some cases; a typical example of this is the Kohn 
mode of oscillation set up by a harmonic trap displace¬ 
ment which is not reproduced by such models m- Vari¬ 
ants of the kinetic model described here, whose single¬ 
component limit does not suffer from such a problem 
[40] . have been put forward in mm-, as explained 
in more detail within the present manuscript, the latter 
work [T161 undertaken by the present authors, was specif¬ 
ically designed in order to introduce the collisional terms 
not explicitly dealt with in previous kinetic approaches, 
in a way which facilitates its numerical implementation. 

The aim of this work is twofold: (i) firstly, we pro¬ 
vide a detailed derivation (Secs. II -IV) of our previously 
proposed multi-component kinetic scheme m, which 
includes both condensate and thermal cloud dynamics 
and all their cross-collisional terms; (ii) moreover, we 
show how numerical application of our scheme to near- 
equilibrium situations (Sec. V) can be used to map out 
regimes of near-hydrodynamic behaviour in accessible ex¬ 
perimental mixtures, clearly highlighting the extent to 
which the relevant degree of ” hydrodynamicity” with re¬ 
spect to different collisional processes can be controlled. 
For completeness, we also briefly describe hydrodynamic 
multi-component equations (Sec. VI) and summarize the 


relevance of our work in the context of existing multi- 
component treatments (Sec. VII). The derivations pre¬ 
sented in the main text are also supplemented by five 
more technical appendices. 


II. COUPLED DYNAMICAL EQUATIONS 


The starting point for our derivation will be the gen¬ 
eral Hamiltonian describing an interacting bosonic binary 
system, with the two components labeled a and b respec¬ 
tively. The Hamiltonian describing the binary system is 
written in second-quantized form as 


H= dr\Y^. 
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and the two-body interactions are given by 

Hr- J dr( Y, 

( 2 ) 

where 4/j = 4/ 7 (r) is the bosonic annihilation operator for 
an atom of species-j, which obey the usual commutation 
relationships for bosons 

[^j( r );^!( r ')] =<M(r-r'), (3) 

['F J -(r),'F fe (r')] = [’J'j(r), ^l(r')] =0. (4) 

The s-wave collisions between atoms in different com¬ 
ponents are encompassed by gkj = 27t hfakj/nikj, where 
akj defines the scattering length between atoms in com¬ 
ponents j and k, and mf} = mf 1 + mZ 1 defines the 
reduced mass. The underlying single-particle Hamilto¬ 
nian appearing in Eq. ([!]) can in general contain exter¬ 
nal potentials, coherent couplings as well as the effective 
trapping and kinetic energies of the atoms. Here Vj( r) 
denotes the trapping potential for atoms of species j and 
can be of any form. 

In the language of symmetry breaking the condensed 
and non-condensed degrees of freedom are separated by 
means of the Beliaev decomposition 


’M r ) = fa( r ) + ^( r )- ( 5 ) 


The condensate of component j is described by the clas¬ 
sical held fj (r) = (’F J (r)), while the non-condensate for 
component j is encapsulated by the fluctuation oper¬ 
ator <5j(r), whose symmetry breaking average is taken 
as zero [117j . i.e. (<5^) = 0. Using the equations of 
motion for the Bose held operators obtained from the 
Heisenberg picture and taking averages with respect to a 
broken-symmetry non-equilibrium ensemble, one obtains 
the equation of motion for component j (for j £ {a, 6}) of 
the condensate held <pj = 4>j(r. t) in the form [ 5811711 . 173 ] 


ih^ = 
dt 
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FIG. 1. (Color online) Schematic representation of the various scattering processes for the binary system. Left: The binary 
system along with all possible transport pathways. Each component is composed of a condensate (below the dashed line) and 
a collection of non-condensate modes (above the dashed line), cumulatively comprising the thermal cloud. Both collisional 
processes (denoted by C and C) and condensate-condensate scattering events that contribute to the mean-field potential, U c , 
seen by each condensate are clearly highlighted. Right: Schematic representation of the coupled equations for the condensates 
(Eq. @) [bottom] and the thermal clouds (Eq. ©) [top] are shown for component a; each diagram represents a momentum 
and energy conserving collision between condensate a ( b ) atoms, shown as closed blue (open red) squares, while thermal a (b) 
atoms are depicted as closed blue (open red) circles. 


+ gkj 
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• R k d=—ig k -{5' h 5j t 8j) / (j)j describes scattering be¬ 
tween different components; 


Here we have defined an effective potential for the com¬ 
ponent j condensate, to encompass, in addition to the 
trap potential, the mean fields of both condensate and 
thermal clouds of both components, via 


• R k ^=—igkj(8l.8j)(f>k/(f>j differs qualitatively from 
the first two (see later) and accounts for an im¬ 
portant ’’condensate collisional exchange” process 
not explicitly included in previous treatments. 


U^(r,t) = Vj{r)+gjj(n c j +2hj) + g kj (n Cyk + n k ), (7) 

where n c j = \<j>j | 2 is the condensate density for compo¬ 
nent j, and hj = hjj = (SjSj) is the (diagonal) non¬ 
condensate density; we also introduce the off-diagonal 
non-condensate density hkj = (<5j<5fc) valid for j ^ k. 
The total density of component j is defined by 


rij = n C}j + hj = \(t>j\ 2 + {Sjdj) . (8) 

Equation can then be written in the simpler form 

dij>j r h 2 




2 in, 


V 2 + U{ - i{l!" + R kj + IT-') 


</>,, (9) 


where the important source terms RR, R k i and R fej ac¬ 
count for atomic transport between the condensate and 
non-condensate for the two components of the gas, and 
are defined in terms of triplet and pair anomalous aver¬ 
ages of the fluctuation operators as: 

• R^=—igjj(8^SjSj}/4>j describes the intra¬ 

component scattering of condensate and non¬ 
condensate atoms, as in the usual single-component 
kinetic equations |3B[, ® , 71, [52] ; 


Within the so called “Popov approximation” (see Ref. 
[71]), pair anomalous terms appearing as ( SjSj ) (diag¬ 
onal) and (8k8j) (off-diagonal) in Eq. § are dropped 
(this is justified on energy conservation considerations - 
see Appendix [B] for a discussion of their physical mean¬ 
ing). 

The corresponding dynamics of the non-condensate 
atoms are encapsulated by coupled quantum Boltzmann 
equations for each component of the gas. Adopting the 
notational shorthand / JJ ( r, p, t) = f J for the distribution 
function, the kinetic equation for component j is written 
as 


|/ j + — p ■ v,,/' - v p /' • v r r ; ; 

ot rrij 

= (^12 + Cl 2 ^ + ®12 + (^22 + ^ 22 ^ • ( 10 ) 


Equation (101 defines the quantum Boltzmann equation 
for component j of the binary system, where each of the 
collision integrals on the right hand side describe qualita¬ 
tively distinct scattering processes occurring within the 
multi-component partially condensed bosonic mixture. 
The non-condensate density associated with component 
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j is given by 


^ (r ’* ) = /(5p /i(r ’ P ’' ) - (U) 


A schematic representation of all arising collisional pro¬ 
cesses for the binary mixture is shown in Figure[l] Equa¬ 
tions © -([T0| represent a closed system of equations, with 
the three source terms of Eq. ([9]) related to the collision 
integrals in Eq. (10) via the relationships 


R kj (r,t) 

R kj (r,t) 
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(12a) 

(12b) 

(12c) 


In this work we first detail the derivation of the above 
equations, which is similar in spirit to the established 
methodology [23 0DJ EH and subsequently use them to 
analyze the relative importance of the emerging colli¬ 
sional processes and the degree of hydrodynamicity of 
typical experimental configurations. 


III. KINETIC FORMALISM 

In order to correctly account for all of the relevant scat¬ 
tering channels amongst atoms in the binary mixture, 
a careful microscopic analysis is required. Pioneering 
work [MUSS] demonstrated how quantum kinetic theory 
could be used to understand the dynamics of the non¬ 
condensate. 


A. Separation of Timescales: Identification of 
Slowly-Varying ’’Master” Variables 

Trapped atoms within the gas are treated as under¬ 
going motion within the trapping potential, which is oc¬ 
casionally interrupted by the s-wave collisions between 
particles. As such, two important collision timescales 
emerge: the duration of a single collision event between 
a pair of particles, which is defined as tq = dkj/v, where v 
is the average velocity of the particles during the collision 
event, and the time in-between collisions t c = 1 /(na^v) 
where n denotes the particle density [68]. As the ki¬ 
netic and interaction energies of the particles are typi¬ 
cally small, the dynamics of the gas are encapsulated by 
a separation of timescales that satisfies tq <C t c , imply¬ 
ing that for weak interactions, we can apply an effec¬ 
tive perturbative treatment, which is fully equivalent to 
the adiabatic elimination of anomalous averages used in 

Refs. fTiminm 

Such an action requires the explicit identification of a 
few slowly-varying ’’master” variables. This task should 
not be taken lightly, as it is the key step determining the 


final form of the equations. By identifying the slowly- 
varying variables, one effectively characterises the mean 
field potentials of relevance in the system (or vice versa), 
thus also fixing the form of the unperturbed Hamiltonian. 
The latter essentially fixes the ” basis” in which the equa¬ 
tions are formulated, i.e. whether one deals with bare 
harmonic oscillator states (as e.g. in [saiMiizs]), dressed 
Hartree-Fock states (as most commonly the case EZ]), or 
even in quasi-particle basis (more challenging, but see 
also 11191. i 120] l. Clearly, all above are inter-related, and 
the importance is to be consistent within a particular 
treatment. In any finite temperature system, we expect 
to have non-negligible components of both the conden¬ 
sate and the thermal cloud of the system: this already 
defines the two slowly-varying quantities as \4>j\ 2 and hj. 

The important question is whether other quantities 
should also be considered as slowly-varying - in this con¬ 
text we should consider the following quantities appear¬ 
ing explicitly in Eq. (|6|: 

i Off-diagonal normal pair averages (<5j<5fc) {j ^ k); 

ii Anomalous pair averages of the form (5jSk) (both for 
j = k and j ^ k); 

iii Anomalous triplet averages of the form (SlSjSj). 

Reflecting from our knowledge of the single-component 
case [67], we note that, as pointed out in EU. the main 
condensate kinetics, i.e. its particle exchange with the 
thermal cloud should come through the latter term. Pair 
anomalous averages could also be included into the treat¬ 
ment through additional self-consistently coupled equa¬ 
tions of motion, as done, for example within the con¬ 
text of a bare particle basis formulation in 68, 69. 73 . 
Their role is discussed in Appendix [B] which shows why 
such terms can, to lowest order, be neglected due to vi¬ 
olating energy conservation. More generally, their inclu¬ 
sion would describe many-body effects mm, which 
are however not expected to be significant in weakly- 
interacting atomic condensates. Based on this, we are 
thus justified in only including such terms in the perturb¬ 
ing Hamiltonian, or even dropping such terms altogether 
from our formalism (the so-called Popov approximation 
[ 21 ). 

This leaves us with the off-diagonal normal pair av¬ 
erages of the form (S^Sk)- In general, these could be 
thought of as describing coherences between the two 
physical systems and could be treated on equal footing 
to condensate and excited state populations pm ms] . 
However, in the absence of any external coupling, one 
would expect such terms to evolve on the more rapid 
collisional timescale, and thus be suitable candidates for 
adiabatic elimination. The fact that they give rise to 
the highly intuitive, but never yet numerically character¬ 
ized, ’’condensate exchange collisional process” confirms 
a posteriori that such treatment was indeed justified. 

Having made such an explicit identification, we can 
now proceed with the perturbative treatment, or adia- 
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batic elimination, of all rapidly-varying off-diagonal nor¬ 
mal and anomalous averages. 

B. Identification of a Perturbing Hamiltonian 

To continue we should now explicitly partition the sys¬ 
tem Hamiltonian given by Eq. 0 as 

H = H M f + (H — Hmf ) = H M f + (13) 

where H' defines the perturbation, and Hmf is the 
quadratic mean-field (unperturbed) Hamiltonian con¬ 


taining only the identified slowly varying quantities (con¬ 
densate mean field and diagonal non-condensate densi¬ 
ties). To proceed, we consider the usual separation of 
the full quartic system Hamiltonian into terms identified 
by a label indicating the number of non-condensate op¬ 
erators appearing in each, i.e. from Hq to Hi (see e.g. 
Refs. |381 mj). This takes the form 

H = H 0 + Hi + H 2 + H 3 + Hi ( 14 ) 

where upon defining h 0 = — (H 2 /2mj)S7 2 + Vj( r) as the 
single-particle contribution from Eq. ([!]), one obtains 


Hn = 


H 1 = 


H-2 = 


* £ 
^ 3 


ho,j + ^-\<Pj \ 2 


fo + ^29kMj\ 2 \(t>k \ 2 k 




% o,j + 




6 j + h.c.^ + E 5 ^' (^*j\(fk\ 2 Si + (fk\(fj\ 2 St + h.c.J |, 
k^j 

s j +^ 2 s}s}+h. c ) +e 9kj A 1 2 5\k + r^JkSi+r^kSis.+h.c. 




H., = 


Hi = 


J dr { Y^9jj (fj&jfijfij + h.c.^ + E_s*tf (tjSWj + tlWk + hx -) }> 

[ A E 9 f¥M +E 9ki¥M 

J ^ 4 U-Ad 


(15a) 
(15b) 
, (15c) 

(15d) 
(15e) 


We wish to work with a reduced unperturbed Hamil¬ 
tonian which is (at most) quadratic, and so we per¬ 
form conventional (but not exact) mean-field approxi¬ 
mations I? 3 mg to only include the leading part from 
the beyond-quadratic Hamiltonian into the unperturbed 


Hamiltonian. Our perturbative treatment of the multi- 
component gas is motivated by Wick’s theorem [123 ]. We 
apply mean-field approximations to H 3 and H 4 above in 
order to reduce these terms to quadratic form. These are 
defined as 


S\.Sk5j ~ + (d\Sj)Sk + (44)<^, 

+<44)^1- [<<5]^)(^4)+(5]4)<^4)+(^^)<44)l, 


(16a) 

(16b) 


valid both for j=k and j^k. The reason we like to work 
with an approximate quadratic Hamiltonian, is because 
this is, at least in principle, diagonalizable by a Bogoli- 
ubov transformation to quasiparticle basis. In what fol¬ 
lows, we do not however consider the dressing of parti¬ 
cles to quasi-particles, but choose to work instead with 
dressed single-particle modes in the Hartree-Fock limit 
E3- Hence, our unperturbed mean-field Hamiltonian 


defining the energy basis of the system ultimately takes 
the form [38U73] ' 

Hmf « (H 0 + 5 H 0 ) + (H 1 +6H 1 ) + (H* ias +6H* iag ). (17) 


Here, the shifts (oc dHfj that appear in each bracket 
are found by applying the mean-field approximations of 
Eqs. ( 16a)-( fl6b| to H3 and H4. The first term in each 
of the brackets in Eq. © describes a contribution from 
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Eq. ([2]) above with the subscript indicating the number 
of fluctuation operators appearing within the operator 
Hi (see Eq. ( 15al-( fl5e| )), while the second term SHj 
arises from mean-field approximations, reducing prod¬ 
ucts of three or more fluctuation operators to quadratic 
form. The ‘diag’ superscript appearing in the final terms 
refer to diagonal contributions with equal component in¬ 
dices. 

The definition of the mean-field Hamiltonian (Eq. 
pT| ) along with Eqs. (15a>-( 15e) and (16al-( l6b| ) then 
allow us to write down the form of the perturbing Hamil¬ 
tonian, a detailed account of which is given in Appendix 

El 


C. Perturbative Description of Condensate and 
Thermal Clouds 


The chosen perturbing Hamiltonian fll(t), (see Ap¬ 
pendix]^ Eqs. (A4a|-(A4d|) will allow us to construct 
our multi-component kinetic theory. It is straightfor¬ 
ward to check that the definitions of H[ (t) along with 
our choice of mean-field Hamiltonian of Eq. (JTt]) recov¬ 
ers the Sclirodinger equation given by Eq. §• 




(18) 


Indeed, it can be seen that the first term on the right 
hand side of Eq. (18) generates the condensate potential, 


mean-field potentials and anomalous pair averages which 
go into the definition of t/|, while the second yields the 
two triplet terms, in agreement with Eq. @. 

To describe the dynamical evolution of the non- 
condensed degrees of freedom, we define the multi- 
component single-particle Wigner operator as mm 


it was assumed that the (matrix valued) non-condensate 
potential U n ( r, t) along with the optical coupling strength 
$~I n (r,t) vary slowly in space, which leads to a qual¬ 
itatively different expression for the kinetic equation. 
Within this approximation the off-diagonal terms in the 
Wigner operator f k i are explicitly computed within the 
perturbing Hamiltonian, leading to matrix valued kinetic 
equations describing the non-condensate dynamics. For 
the two cases of optically coupled condensates with either 
spin-1 or spin-1 internal degrees of freedom, the relevant 
kinetic equations are given by Eqs. (52) and (41) in Ref. 
m and (1151 respectively. We are however interested 
in understanding an incoherent binary mixture, hence 
for j ^ k we set in the final calculations / fcj ( r, p,t) = 0, 
i.e. no explicit long-lived coherences between off-diagonal 
normal pair averages. The Wigner operator directly al¬ 
lows us to calculate relevant non-equilibrium expectation 
values for the multi-component system. The correspond¬ 
ing equation of motion for the diagonal elements of the 
phase-space distribution function / J (r,p,t) is written 

^f (r, P, t) Tr p(t, t 0 )[f (r, p, t 0 ), #mf(£)] 

+ (r,P,to), H' (t)\. (22) 

where the first term on the right hand side gives the 
free streaming terms and the second term generates the 
individual collisional integrals. 


IV. DERIVATION OF COLLISIONAL 
INTEGRALS 


A. Mathematical Formalism 


f kj = j dr'e ip r '/^t( r + r'/2,t)6 k (r - r'/2 ,t), (19) 

where the corresponding phase-space distribution func¬ 
tion is defined as (r, p, t) = Tr p(t, to)/ fej (r, p, t 0 ), and 
pit, to) defines the general density matrix of the system, 
which is related to the initial density matrix p(to) by the 
unitary transformation pit,to) = U(t, to)p{to)U^(t, to). 
The unitary evolution operator satisfies the equation of 
motion 


intuit, to) = H(t)U(t, to) (20) 

ot 

while the density matrix p(£,£ 0 ) evolves according to 

d 

ihg- t Pit, to) = [Hit ), pit , to)]- (21) 

In general the coherences of the non-condensate are non¬ 
zero only when an optical or magnetic coupling exists be¬ 
tween the states |a) and |5). This particular case was ex¬ 
plored in 1141 1T51 for spinor Bose gases. In those works 


In order to calculate a closed set of equations describ¬ 
ing the finite temperature dynamics of the bosonic mix¬ 
ture, the collision integrals appearing in the dissipative 
Schrodinger equation (]9j) and the quantum Boltzmann 
equation (101 are derived using the perturbation Hamil¬ 
tonian, H' , defined by Eqs. (A4a|-( A4d) in Appendix [a] 
This is in turn accomplished by expanding the fluctua¬ 
tion operators in terms of their Fourier components, and 
calculating the non-equilibrium expectation values of the 
various products of such operators. The non-equilibrium 
average of an arbitrary time-dependent operator O(t) can 
be computed using the general density matrix p(t, to) de¬ 
fined previously, as well as the mean-field evolution op¬ 
erator Sit, to) that satisfies the equation of motion, 


ih^Sit, t 0 ) = H M Fit)S(t , t 0 )- (23) 

It can then be shown that the expectation value of the 
operator O(t) can be written as [57] 


{Ot) — Trp to < S tM O t0 St, tl 
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dt'sl l:to [sl t ,d to s t ,t’, 


H't,]Si 


t'yt 0 


to 


(24) 


where A tlit2 = Aftijt?) has been used here and in what 
follows to abbreviate the time dependence of time evolu¬ 
tion operators. We use Eq. (241 to compute closed ex¬ 
pressions for the source terms appearing in Eq. (|9[ ). along 
with the collision integrals appearing in Eq. (|10| ) above. 
The first term on the right hand side of Eq. ( |24| can be 
dropped, as it is assumed that for long times any initial 
correlations present in the system vanish (the Markov 
approximation). Thus, Eq. (24) becomes 


tors, which for component j is given by 


^■m°) = 


*3, P c 


Jp-r/h 


(29) 


The expansion defined by Eq. ( |29[ ) allows us to write 
the Fourier transform of the Wigner operator defined by 
Eq. (19) above. This is best handled by switching to 


the center of mass and relative momenta for the two- 
component system. Hence, the general Wigner operator 
in momentum space for a binary mixture is written as 


/^(r,p,to)=e i 2 l ^ r 'P/ s 




y,2(mj /M)p-q/2 a fc,2(m fc /M)p+q/2 e 


-q /n 


(30) 


(6 t ) ~ -- / dt\Sl^S{ t ,6 t0 S t ,t-,H[,]S tl , t0 ), (25) 


*0 


where (...) = Tr/5 to (...) for the right-hand side. 

Since we have identified the condensate and non¬ 
condensate fields as slowly varying, we write n c j(r', t 1 ) — 
n c ,j(r,t),hj(r',t') ~ fij(r,i) and Ul(r',t') ~ U£(r,t). It 
is useful to write the condensate wave function for com¬ 
ponent j in the density-phase representation using the 
Madelung transformation <f)j = exp(i9j), in which 

case can be expressed as 


and the total mass is M = rrij+mk■ Since we are only in¬ 
terested in the (incoherent) processes involving the diago¬ 
nal elements of the Wigner operator, (see Refs. [11411115 
for generalizations that include the off-diagonal contribu¬ 
tions to the Wigner operator.) we work in what follows 
with Eq. (30) in the limit j = k. Hence 


/ J (r,p,t 0 ) = XXp- 


q/2 a hP+q/2 e 


-q/fi 


(31) 


Equation (31) will be used to calculate the collision inte¬ 


grals in the following sections. 


0j( cs 6j(r, t) + -^(t>-t) + V9j ■ (r' - r), (26) 


- @j( r i t)~ 




r'-r). (27) 


In writing Eq. (27) we have used the Euler equation 


for component j (see Eq. (80) in Sec. VI) in order to 
introduce the local condensate energy, 


1 -2 

£ J c {r,t) = Mc( r ,t) + -mv> c . 


(28) 


Finally, the Fourier transform of H'(t) allows us to derive 
non-equilibrium expectation values for arbitrary prod¬ 
ucts of operators in the following subsections. In order 
to calculate closed expressions for the collisional inte¬ 
grals appearing in Eq. © and (fl0|), we must express 


the higher order correlation functions (those formed from 
non-equilibrium expectation values of products of fluctu¬ 
ation operators) in terms of the distribution functions 
/ J (r, p,f). As such, the perturbing Hamiltonian H'it) 
is used to extract collision integrals to second order in 
the scattering length akj, while maintaining the effect of 
interactions in the collective mode energies and chemi¬ 
cal potentials to first order in a^j . in the spirit of the 
single-component ZNG approach [124] . 

The evaluation of non-equilibrium quantities requires 
the Fourier expansion of the non-condensate field opera- 


13. Source Terms from Anomalous Averages 

1. Condensate Growth Terms R 33 and Rr 3 
(from triplet anomalous correlations) 

We begin by considering the triplet contributions to 
Eq. (|9|, R V and R k3 . We will explicitly compute R kl , 
using Eq. (24). The triplet contributions can be decom¬ 
posed as 

(SiSkSj) = (S f M)(0 A i^k^kSj)( 3 ), (32) 


The two terms in Eq. (]32j) above require the computation 
of averages from the perturbation Hamiltonian involv¬ 
ing on e [Eq s. (A 6 a)-(A 6 bl) for j=k] and those for three 
[Eqs. (A 6 el-( A 6 f| ) for j^k\ fluctuation operators respec¬ 
tively. We first compute (SjSjSk)( 3 ), be. those contri¬ 
butions arising exclusively from commutations involving 
the perturbing Hamiltonian As such we first 

calculate 


( 3 ) = 


~n J dt '& 

to 


[SlJlhSjSt, 


t',t 0 L u t,t' u k 


! H'a,kj]St',to) ■ 


(33) 


After using the definition of the Fourier transform of 
H'S’kjit) defined as 































Kkj 


«> = ^£ E 

k^j P2.P3.P4 


^Pc + P2,P3+P4 e 


-t)/n-p j 0 -r/h) - i 


k, p 2 afc .P3 a j,P4 


+ h.c. 


(34) 


r 


Eq. (33) becomes 


<Wi> ( 3)=-|^i E *(4+4-4-4) 


X(5 ; 


P2,P3,P4 
fkf fk | i W yd 


Pi+P2,P3+P4 


[/£(/£ + !)(/! + l)-(/ 2 fe + l)/ 3 Vi], (35) 


where the shorthand f k = / fc (r, p„,f) has been used 
in the above and what follows, e J c = p J c + ^rrijV^ j and 


£ p = P 2 /^ m j + 4 define the non-equilibrium condensate 
and thermal energy for atoms in component j respec¬ 
tively. By writing Eq. ( |35| we let to —> oc in order to 
evaluate the integral over t in Eq. 


^ (See Eq. |C19 ) 
and discussion in Appendix [C] for an explanation of this 
important step). Calculation of (SlSkSj}( i) allows us to 
write the first term in Eq. ( [33] ). Then by taking the con¬ 
tinuum limit, we obtain the source terms 


R kj = 


R n = 


4 


2(2tt) 5 H 6 

9 2 j 

( 2tt) 5 Ii 6 


dP2 / dP3 / dP4 <^(P^+P2-P3-P4 )${ £3 c+ £ p 2 - £ p 3 - £ L) 


dp2 / dp3 / dP4 S( pJ+p 2 -p3-P4)(5(e^+eL-ep3 _e P 4 ) 


/2 fc (/3 fe + l)(/I+l)-(/2 fc + l)/ 3 fc /l 


fiui+wi+v-mmfi 


(36) 

(37) 


Here, Eq. (37) has been obtained by repeating the steps 
following Eq. (pi) for R k ' J . 


pair of fluctuation operators, {S\Sj). Hence we calculate 


2. Condensate Exchange Terms Et fe -' 
(from off-diagonal normal pair averages) 


< 3 &>=-i / ( 38 ) 


to 


The final dissipative source term IR, fcj appearing in Eq. 
is comprised of a normal average of an off-diagonal 


Computation of Eq. (38) requires the Fourier transform 


of H’ 2 k -(t), the relevant contribution being 


^2,yW E E 9kjy/ n c,j n c,k^ *pi+pj,P 2 +P* e 

kjtj P1.P2 


-i« e ,- ek )-(4-eW-t)/H-(P°^y^al pi a^ P2 + h.c. j>, (39) 


r 


Using Eqs.(38) and (39) yields the expression 


= -ygwMt E *(4 + 4! - 4 - 4) 

Pl,P2 

X 4 +Pll pj +P2 [ifl + l)/2 - flirt + !)]• (40) 


By taking the continuum limit of Eq. (40), the off- 
diagonal pair average becomes 


T3 kj _ 9kj 

2(2tt y-h? 


n c ,k / dp I / dp-2 S(p k c +p 1 -p J c -p 2 )S(£ k + £ J P1 -£ J c -£p 2 ) 


ifi + m- flirt + i) 


(41) 


The three expressions derived in this section, Eqs. 
(36), (37) and (41) are the important source terms that 


appear in the dissipative Schrodinger equation. Equa- 
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tion (41) arises due to our explicit separation of slowly by Eq. (34). The multi-component nature of the prob¬ 


and rapidly varying quantities in the system Hamilto¬ 
nian, and can be understood as a collisional energy ex¬ 
change process between the two condensates, whereby 
a condensate and thermal atom in differing components 
scatter into corresponding thermal and condensed states 
respectively. 


C. Quantum Boltzmann contributions 


1. Collisional an d Cy| terms 


To complete the derivation, we require the collision in¬ 
tegrals appearing on the right hand side of Eq. ©• 
These are computed using the definition of the multi- 
component single-particle Wigner operator (Eq. ©), 
along with the Fourier transform of H' 3 fc -(i), as defined 


lem leads us to partition the “C 12 ” collision integral into 
two parts, the first Cy 2 defines the intra-component scat¬ 
tering of atoms, while the second Cyj gives the inter¬ 
component collision rate. We wish to calculate both 

C 12 = ^Trp(t,to)[p(r,p,to),H' 3J (t)], (42) 


and 


C i 2 = Ptt’toW {r,P,t 0 ), H' 3kj (t)]. 


ih 


(43) 


As before, let us illustrate the derivation of these 
two terms by computing the off-diagonal contribution, 
Eq. (43). This is accomplished by using Eqs.(34) and 


(43), giving 


r< k i — 
°12 — 


9kj 


ihW 




; PJ+P2,P3+P4 < ' 


&k ( ^P2,P+q/2(®j !p _q/2%,P3®fe,P4) <^p 3 .p-q/2 


x a J 


,p+q/2«fe,pi)^ Vc+P 2 ,P 3 +P 4 43 ^P4,p-q/2 <ai,p 2 “fc,P3 ®j,P+q/2) ll.C.|. 


(44) 


r 


Then, by using the definition of the multi-component 
three-field correlation function, the continuum limit can 


be obtained as before by replacing the summations with 
integrations, giving 


c k i — 
°12 — 


9li 


(2t r) 2 fi 4 


n c ,k / dp 2 / dp 3 / dp 4 S(p k + p 2 - p 3 - P4)<5(£c + e p 2 “ e P3 “ £ P4 ) 


(ti + vtifz-mi + mz + i) 


<5(P - P 2 ) - <5(p - P 3 ) 


J dp 2 J dp 3 J d Pi S(pi + p 2 - p 3 - P4 )5(e J c + £ p2 - e p3 - e 3 Pi ) 


(f k +1 )fZfi - / 2 fc (/ 3 + i)(/i +1) 


(~<33 _ ^9jj 
12 "(2 tt)W U c ’° 


S{p~Pi), 

rim I dp-2 I dp 3 I dp 4 S(p J c + p 2 - p 3 - P4 )S(e J c + £ J P2 - e J P3 - e° Pi ) 

<5(P - P 2 ) - <5(p - p 3 ) - d(p - p 4 ) 




(45) 


(46) 


with Eq. (|46j) obtained by repeating the same steps for 
the collisional integral defined by Eq. (42). It can be seen 
that Eq. (46) is equivalent to the C 12 collision integral 
from the single component kinetic theory mm 


2. Exchange collisional term 

To complete our discussion of collisions involving con¬ 

densate and non-condensate particles, we compute the 
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exchange collisional integral, which is defined by 


us to compute an expression for given by 


C i2 = ^Trp(i,t 0 )[/ J (r,p,to),-ff2,fcjW]- (47) 


_ 2 9kj 

12 ih 


E Vc+Pl,P^+P2 ( ^'^fc( a ip 2 Sfe .Pl) 


Pi 5P2 




P2/ 


(48) 


Following the methodology discussed above, the Fourier 
transformed Wigner operator along with Eq. (39) allows 


upon inserting the pair correlation function into Eq. (48) 
and taking the continuum limit yields the exchange inte¬ 
gral C& given by 


,kj _ 2n 9lj 


C 12 = n c ,k n c j / dpi / dp 2 d(p^ + pi-Pc-p2)<5(^ + ep 1 -£c- e p 2 ) 


/fC/2 + l)-(/f + l)/i? S(p-p 2 ). (49) 


3. Collisional Coo and Cll terms 


The final collisional processes described by Eq. (10) 


are the terms that account for interactions exclusively 
between non-condensate atoms, C 2 i and C 22 ■ The quan¬ 
tities we wish to evaluate are given by 

C 22 = ^Trp(i,t 0 )[/' , (r,p,to),H'lj(i)] ! (50) 


and 


^22 = P&toW {r,P,t 0 ), Hi, kj (t)]. (51) 


iH 


As with the Cfl and C ^ collision integrals, the commu¬ 
tation of the Wigner operator p (p, r, t) with H' 4 j(t) and 

H 4 k j{t) generates the intra ( C%£) and inter (C 2 l) colli¬ 
sion integrals respectively. We illustrate the derivation 
by considering the collisional integral C 22 ■ This requires 
the Fourier transform of H 4 fej (t), which is given by 


-®4 ,kj 9kj 

ko K Pl.P2.P3.P4 


XX ^Pl+P2,P3+P4®j |Pl ®i,p 2 ®fc,P3®j,P4 ^jk 


E«px. At 


P2 a j,pi a fe,P2 


P1.P2 


71 kj XX ^Pi'P2^fc,pi®i.F 

P1.P2 


(52) 


r 


By inserting Eq. (52) into Eq.(51) and taking the contin- 

~<kj 
'22 


uum limit, the collisional integral C 2 % is found to be 


fykj _ 9kj 
22 (2n) 5 h 7 , 


f<n_ 


dpi I dp 3 I dp i d{p+p 2 -p 3 -pP)5{£ : ‘+e k -£ k -£ 3 ) 


(27 t) 5 H 7 , 


dp 2 / dp 3 / dp i 8{p+p 2 -p 3 ~pp6{e J +E 7 -e J -el ) 


(/ j + l)(/2 fc + l)/3 fc /i-/V 2 fe (/3 fc + l)(/I + l) 

(/ J '+i)(/l+i)/l/i-/V^(/l+i)(/i+i) 


(53) 

(54) 


Equation (54) is obtained by repeating the above steps 
for the collisional integral defined by Eq. (501. This for¬ 


mally completes the derivation of all collisional integrals 
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appearing on the right hand side of Eq. (101. The expres¬ 
sions given by Eqs. (45), (46), (49), (53) and (54) will be 
used in the subsequent sections to study the equilibrium 
properties of binary condensates. 


V. NUMERICAL SIMULATIONS 


At finite temperatures, the various collisional processes 
have a heavy influence on the coupled dynamics between 
the condensates and the non-condensed atoms, such as 
the damping of collective modes |125H127j . the decay of 
solitons m and vortices m , as well as the growth 
of the condensate [82, [83], as seen in the correspond¬ 
ing single-component kinetic theory. In this section, af¬ 
ter obtaining our equilibrium distribution (Sec. VA), we 


compare the roles of different collisional processes under 
the variation of isotropic trap frequencies (Sec. |VD 1 and 
trap geometries (Sec. VE). T his is achieved by calcu¬ 
lating the collisional rates (Sec. VB) and hydrodynamic 
parameters (Sec. V C) for various equilibrium binary sys¬ 
tems. We show and explain the scaling relations be¬ 
tween the hydrodynamic parameters and isotropic trap 
frequency in Sec. |VD| In Sec. |VE[ we demonstrate the 
generic dominance of the exchange collisional process C 12 
across the different trap geometries, even though all col¬ 
lision rates strongly depend on the relevant scattering 
lengths. Importantly, our results in Sec. |VD| and |VE| 
illustrate different ways to control the hydrodynamicity 
of the collisional processes, which can be of high interest 
to experiments. These include 


or 51ao [TB] 1 enable the probing of both miscible (A = 
< 712/^511522 < 1) and immiscible (A > 1) regimes. In 
a previous work m, we have presented our numerical 
results for such systems in an isotropic harmonic trap 
(frequency u> = w_l = ui z = 2n x 20Hz). In particular, 
we have highlighted the dominance of the exchange col¬ 
lisions C 12 over the other collisional processes within the 
temperature range 0.3 < T/T c < 0.9. Here, we perform 
a more detailed analysis that compares rates for different 
isotropic trap frequencies and different trap geometries, 
using our results for uj = 2tt x 20Hz as a reference. 


A. Equilibrium solutions and condensate fractions 


The equilibrium density distributions at temperature 
T are numerically obtained by setting the source terms 
(R 33 , R k3 , ) and the collision integrals (C 12 , C 12 , C 22 ) 

to zero and self-consistently solving Eqs. ([£]) and ( fTo| . In 
order to speed up the computation, we adopt the semi- 
classical approximation m for the local non-condensate 
density 


nj{T,t)=j = ^g s/2 ( Zj ), (56) 

where A j = ^2nh 2 / (rrijkBT) is the thermal de Broglie 
wavelength, Zj( r) = exp{[fj, 3 c — Ul(r)]/(kgT)} is the local 
fugacity and the chemical potential n 3 c is obtained from 
the imaginary-time evolution of the condensate equa¬ 
tion Q. 


(i) bringing the hydrodynamic parameters of the vari¬ 
ous processes closer in magnitude by increasing the 
trap frequency and the temperature, 

(ii) increasing the hydrodynamicity of all processes to¬ 
wards the hydrodynamic regime by changing trap 
geometry, 

(iii) controlling the hydrodynamicity of the intraspecies 
and interspecies collisional processes by tuning the 
relevant scattering lengths through inter- or intra¬ 
species Feshbach resonances. 

In the final subsection |VF[ we briefly explore the validity 
of the usual high-temperature approximation (/?(e — /j) -C 
1) in the context of collisional rates, specifically for the 
exchange collision € 12 . 

Our numerical analysis focuses on experimentally rele¬ 
vant equilibrium 87 Rb- 41 K and 87 Rb- 85 Rb mixtures with 
a total atom number Nj = 10 5 in each component 
trapped in harmonic potentials 

Vj(r) = ^ (u 2 ± (x 2 +y 2 )+u 2 z 2 ). (55) 

These mixtures were chosen as their tunable scattering 
lengths (aRb87 = 99ao, aK4i = 60ao, <XRb87-K4i = 20ao 
or 163ao HU H29]; a Rb87-Rb85 = 213ao, ORb85 = 900ao 


h 2 
2 rrij 


v 2 + u° c 


/4<Ar 


(57) 


We start our analysis by first considering the conden¬ 
sate fractions of the binary mixture at different tempera¬ 
tures T, as shown in Fig. [2] While our method is strictly 
not valid for T close to the critical temperature T c due 
to critical fluctuations, we can nevertheless extract T c 
by fitting the fractions with N c /N = 1 — (T/T c )“ [ 1321 
and compare the extracted T c to the expected shift in T c 
due to finite-size corrections [ 133) and mean-field correc¬ 
tions [134] . 

For a single-component Bose gas and using our 
simulation parameters, T c decreases by approximately 
0.73 AU 1 ^ 3 =2% due to the finite number of atoms. The 

mean-field shift, calculated by —1.3(a/aho)A/ ( 1 ^ 6 , where 
a is the relevant scattering length and aho is the rele¬ 
vant harmonic length, further decreases our T c by 1-2%. 
However, for the T c of 85 Rb in the miscible mixture, the 
mean-field shift amounts to approximately 17% due to 
the large scattering length aRb 85 -Rb 85 = 900ao- Note 
that we did not take into account many-body effects be¬ 
yond mean-field theory 135 , which can instead increase 
the critical temperature. Overall, our extracted T c are 
close to the mean-field predictions for a single-component 
gas [HI]. 
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FIG. 2. (Color online) Condensate fractions of 87 Rb- 41 K (top) and 87 Rb- 85 Rb (bottom) mixtures at different temperatures 
in an isotropic harmonic trap (trap frequency oj = 2n x 20Hz) with scattering lengths aRb87-Rb87 = 99ao , aK 4 i-K 4 i = 60ao 
, dRb87-K4i = 20ao (miscible) or 163ao (immiscible O 129], dRb87-Rb85 = 213ao, and dRbss-Rbss = 900a 0 (miscible) or 
51ao (immiscible) iT^; each species has a total of IV = 10 5 atoms. Dashed lines give the prediction for non-interacting single¬ 
component trapped gas, with condensate fraction N c /N = 1 — (T/T 3 ) 3 , with critical temperature = 42nK. The solid lines 
is our numerical fit using the condensate fraction N c /N = 1 — (T/T c )“, where the extracted T c (indicated in the legend) are 
lower than the mean-held single-component T c [m| by at most 5%. 


87 Rb- 41 K 



7'/ !}io 


FIG. 3. (Color online) Miscible (left) and immiscible (right) 8 'Rb- 41 K mixtures in an isotropic harmonic trap (frequency 
uj = 2tt x 20Hz) at temperature 21nK. Other simulation parameters (scattering lengths and total numbers of atoms) are the 
same as Fig. 2 The reference harmonic length is £ho = \/fi/(m.Rb 87 d;). Top: Condensate and thermal densities. Middle: 
Spatially resolved collision rates between condensate and thermal atoms. Bottom: Spatially resolved collision rates between 
thermal atoms. 


Typical density profiles of binary systems are shown 
in the top panels of Figs. I ( 87 Rb- 41 K) and g ( 87 Rb- 
85 Rb), where the two condensates (dashed lines) mix 
(left columns) or phase-separate (right columns), but al¬ 


ways sit on top of more diffused non-condensate clouds 
(solid lines). These non-condensate clouds have long tails 
that extend much further than the condensate clouds, a 
feature that is also seen in the single-component Bose 
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87 Rb- 85 Rb 


Miscible (A = 0.7) Immiscible (A = 3) 



r/ 4o 


r/ 4o 


FIG. 4. (Color online) Same as Fig. [ 3 ] but computed for 87 Rb- 85 Rb mixtures. 


gas 1271 132j . However, in contrast to the single-peak 
structure in a single-component gas, mean-field repul¬ 
sion from both condensates in a binary mixture leads 
to a double-peaked thermal structure at the condensate 
edges, where the effective mean-field potentials of the 
non-condensed atoms are local minima. The middle and 
bottom panels give the spatial collisional rates involv¬ 
ing collisions between thermal-condensate atoms (II 2 
and Tp) and thermal-thermal atoms (T 2 2 ), respectively. 
These spatial rates depend strongly on the condensate 
and thermal cloud density profiles. In the next subsec¬ 
tion, we give more details on the calculation and analysis 
of these collisional rates. 


B. Collisional rates 


With the equilibrium density profiles, we can proceed 
to evaluate the local collision integrals (45), (|46|) , (49), (53) 


and (54) (where the condensate energy e° c = p J c and the 
condensate momentum p J c = 0 at equilibrium). Since 
these integrals are identically zero at equilibrium, we re¬ 
express them in the form 


C kj out 

12 — °12 


-c 


kj, in 


(58) 


(and analogously for C 12 and C 22 ) to explicitly identify 
the “in” and “out” scattering rates. In this way, we can 
assess the importance of the various collisional processes 
by comparing either their “in” rates or their “out” rates. 
By integrating over momentum space, we obtain the col¬ 
lisional rate 

T^fc/.out f dp si hi. out / - n \ 

r !2(22) = J ^ 12 ( 22 ) (59) 

that measures the number of non-condensed atoms that 
have collided through a particular “out” process per unit 


volume per unit time. The mathematical steps needed 
to compute Eq. (59) are given in Appendix [D] In the 
following, we only give the final formulae used in our 
numerical computation. 

To calculate the collision rates between non-condensed 
atoms (for both k = j and k 7 ^ j), it is convenient to 
transform to the center-of-mass frame. We therefore ob¬ 
tain 


Y'kj, out [ dp 1 tj f dp 2 tk 

22 J (27 Thf h J {2nhf h 

x / V2 K/s fe + !)(/i + !). (60) 


where a k j = (1 + Skj^n a)h is the cross section, Vi and 
v 2 are the initial velocities of atoms j and k respectively, 
H specifies the solid angle of the final relative velocity 
v 4 - v 3 . 

For collisions between condensate and non-condensate 
atoms, we first look at the C-f| process (for both k = j 
and k 7 ^ j) which is present even in a single-component 
Bose gas. We specifically evaluate the “out” collision rate 
that represents the scattering of a non-condensed atom 
from a condensate to produce two non-condensed atoms, 


y\kj, OUt 
1 12 



d P2 fk 

( 2 ?rh ) 3/2 


n c,j 


G kjV 


out 

r 



where w“ ut = |v CJ - — v 2 | 2 — 2{Ui — H 3 c )/mkj is the rel¬ 
ative speed of the initial states corrected to take into ac¬ 
count the local conservation of energy. The reverse pro¬ 
cess, where two non-condensed atoms collide such that 
one of them goes into a condensate, is given by the “in” 
rate as 


jU-yin 
1 12 


dp4 n Ctj <T kj (m k /m k j) 3 
(2nH) 3 4 47 t (1 + mj/m k )\v'™\ 


dv/ 3 fc , (62) 
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where v™ = V 4 — v J c is the velocity of thermal atom j 
relative to the local condensate velocity while the second 
integral is a two-dimensional integral over the velocity 
vector v which is in a plane normal to v“. The velocity 
of the other incoming thermal atom v| is then given by 


fc ■ in (Uj-ni)*? 

3 2 r + + mj |vj?| 

and the outgoing velocity of the thermal atom is 

by 


(63) 

given 


v 2 fe = K/m fe K n + v 3 fc . (64) 


Note that we follow 1271 and drop the cubic term 72 / 3/4 
in numerical simulations as it cancels exactly between the 
“in” and “out” rates. 

Finally, we consider the exchange collisions (M 
j only) novel to our treatment of the binary Bose gas, 
which describes a process whereby one condensate atom 
(say atom k) collides with a non-condensed atom j and 
are then scattered into a thermal (atom k) and condensed 
(atom j) state. The collision rate is 


pout _ 

1 C — 


-M-kj 
771 kj 


77c : k 7lc,j7J r 


47r 


mi+ 1 ), ( 65 ) 


where = m k 1 — rrij 1 plays the role of an effective 
reduced mass while the effective relative speed is 


v r = 


\/Wc,j- v c , fc | 2 -2([t/*-^]-[t/i(-^])/A4 fej . (66) 


Equations (60), (611 and (|65j) are the key quantities com¬ 


puted in our simulations once the equilibrium densities 
are obtained. The first two are evaluated using Monte 
Carlo sampling of the integrals while an exact expression 
(see Sect. JVF| for detailed discussions) can be obtained 
for Eq. (65) because v C J = v C}k = 0 at equilibrium. 

Examples of these spatially-resolved collision rates are 
shown in the middle (Ti 2 and T®) and bottom (T 22 ) pan¬ 
els of Figs. [3] and [4] Note that F 12 are drawn on a larger 
scale compared to r 22 . Typically, the rates between con¬ 
densed and non-condensed atoms (both Ti 2 and T®) fea¬ 
ture localized peaks at the condensate edges where the 
thermal cloud and condensate overlap the most, while the 
rates between non-condensed atoms (r 22 ) follow closely 
the shape of the thermal density profiles. This is because 
r 12 and Tc are approximately proportional to the prod¬ 
uct of condensate densities and thermal cloud densities 
while r 22 is approximately proportional to the product 
of two thermal cloud densities. These observations are 
important to understand the variation of these collision 
rates with respect to the trap frequency; see Sec. |VD| 

On the other hand, comparison among the C 22 or Ci 2 
processes shows that the relative peak values of F can 
be estimated by the relevant cross section a oc a 2 s . For 
example, s 'Rb intraspecies collisions (red dashes) and 
the 8 'Rb- 41 K interspecies collisions (black dots and black 


dash-dots) have comparable peak heights in the immisci¬ 
ble case (Fig. [3]) because of similar cross sections, whereas 
the 85 Rb intraspecies collisions (blue thick dashes) dom¬ 
inates over both 8 'Rb- 85 Rb interspecies collisions (black 
dots and black dash-dots) and 87 Rb intraspecies colli¬ 
sions (red dashes) in the miscible case (Fig. [4]) because 
of the large 85 Rb scattering length a Rb85 = 900a o . 

Finally, for the case of a 8 'Rb- 85 Rb mixture, the sharp 
peak for the interspecies exchange collision Fc (green 
curves in the right middle panels of Fig. [4]) is a con¬ 
sequence of the small mass difference between the two 
different atomic species. This is most easily seen if we 
consider a small spatial region around a critical radius 
r c , at which v r = 0. In this case, we can approximate 
v r = \/C(r — r c )/ A4fc7 f° r some constant C and substi¬ 
tute this into Eq. ( [65] ) to show that the spatial width 
of Tc is proportional to . For small mass differ¬ 

ence, while Tc is sharply-peaked in space, it is neverthe¬ 
less possible to obtain the total number of interspecies 
exchange collisions, a physically meaningful and exper¬ 
imentally relevant quantity, by integrating F® over the 
full cloud volume. 


C. Hydrodynamic analysis 


From the collision rates, we can further extract the 
mean free time r as [3D] 


1 

r 


JVcoll 


dr r(r), 


(67) 


where jV co h is the relevant number of available non- 
condensed atoms taking part in collisions for each pro¬ 
cess. For example, with F^ in Eq. (61), 


N coii = J = J dT 


( 68 ) 


In the case of thermal-thermal collisions (C 22 processes) 
that involve two different components, N co n refers to the 
number of non-condensed 8 'Rb atoms. This choice has 
no significant impact as, for our simulation parameters, 
the condensate fraction of 87 Rb differs from the fraction 
of the other component ( 41 K or 85 Rb) by less than 10%, 
except for the miscible 87 Rb- 85 Rb mixture due to the 
strong mean-field corrections to the 85 Rb fraction, see 
Fig. [2] 

We would like to mention that one could also use an 
alternative time scale defined by t [1241 

i = [dr F(r), (69) 

T IVcon J 

where N con is the relevant number of condensed atoms. 
The key difference lies in _/V C oii and N con . which simply re¬ 
flects our interest with respect to the change in the num¬ 
ber of either non-condensed or condensed atoms. These 
two time scales therefore differ by the order of condensate 
fraction. 
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FIG. 5. (Color online) Variation of hydrodynamic parameters r] = 1/(wt) of 8 'Rb- 41 K (left) and 87 Rb- 85 Rb (right) mixtures 
with respect to the scaling of trap frequency and temperature, (w, T) = k(uio, To) for C 22 (top), C 12 (bottom) and C 12 (bottom, 
green inverted triangles) processes. Scattering lengths and atom numbers are the same as Fig. [5] and [I] The reference frequency 
and temperature are wo = 2-7T x 20Hz and To = 25nK. Each data set of ?; has been normalized by its value 770 at k = 1. The 
solid lines give our predictions ( |72| ) for a = 1 (top) and a = 1/2 (bottom). l/(wrc) departs from our prediction (a = 0) in the 
miscible 87 Rb- 41 K case (first column bottom) because our assumption of localised C 12 process is no longer valid; see Fig.[3]for 
an example. 


When r is compared with the trap frequency, which 
governs the oscillation frequency of a collisionless classi¬ 
cal particle in the harmonic trap, we obtain the dimen¬ 
sionless hydrodynamic parameter 



(70) 


If r/ > 1, a non-condensed atom will experience on aver¬ 
age at least one collision before completing an oscillation 
in the harmonic trap, hence the system is in a hydrody¬ 
namic regime. Otherwise, the system is in the collision¬ 
less regime. These hydrodynamic parameters are highly 
relevant to cold-atom experiments as they determine the 
thermalisation rates [136H139] . In particular, the inter¬ 
species hydrodynamic parameters are crucial to the effi¬ 
ciency of sympathetic cooling [137] . Understanding these 
parameters can therefore help to optimise future studies 
of the various cooling stages. 

In the following subsections, we analyse the collisional 
processes using the hydrodynamic parameter given by 
Eq. ([67). 


D. Trap frequency variation 

Our previous work m has shown that the hydrody¬ 
namic parameter 77 c = l/(wrc) of the exchange process 
can be one to two orders of magnitude larger than the 
corresponding parameters of the C 12 and C 22 processes. 


In this subsection, we show that by varying the isotropic 
trap frequency, it is possible to bring 77 of the various 
processes closer in magnitude. We provide a further ex¬ 
planation based on the scaling of length, energy and con¬ 
densate densities. 

In order to make meaningful comparisons, we scale the 
trap frequency w and the temperature T simultaneously 
by the same factor n such that the condensate fractions 
of the binary mixtures remain approximately the same 
as oj and T are varied. This can be easily understood for 
the non-interacting Bose gas, where the single-particle 
energies appear as multiples of huj and the thermal oc¬ 
cupation (determined by the ratio fno/ksT) thus remain 
unchanged. 

We use Wo = 27 t x 20Hz and T 0 = 25nK as references 
for trap frequency and temperature and consider four dif¬ 
ferent sets of isotropic trap frequency and temperature, 
(w,T) = k x (wo,T 0 ),fc € {1,2,4,8}. The numerically- 
obtained hydrodynamic parameters are shown in Fig. [5] 
which clearly shows that 

- = (71) 

Vo 

where 

r 1 (C22) 

a « < 1/2 (C 12 ) (72) 

{ 0 (C 12 ), 

and 770 = l/(wo t 0 ) is the reference hydrodynamic param¬ 
eter (different numerical values for different processes) 
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while To is the mean free time of the various processes 
at k = 1. In order to understand Eq. (72), we rewrite 
Eqs. (60), (61) and (65) in terms of dimensionless vari¬ 
ables for position r = r /£, momentum p = p t/h and 


velocity v = \/(£oo), where t = \Jh/{rnnbOo) is the har¬ 
monic length for 8 ' Rb atoms and we choose the mass of 
87 Rb here simply as a reference. For the C 22 processes, 
we have 


1 1 f dr dpi dpi 

= pnzzJ ( 2 tt )6 

x J - v 2 | flfzifi + 1 )(/I + !)• (73) 


Since the phase-space distribution / as a function of r 
and p remains approximately unchanged as we vary /c, 
the sole dependence on k in l/(wr 2 2 ) comes from the pref¬ 
actor l/£ 2 oc u> oc k and thus a = 1 , a result consistent 
with those from m- 

We perform the same transformations on the C 12 col- 
lisional process and arrive at 


1 _ 1 

wr i2 ^con 


dr dp 2 fk 
(2tt) 3 h 


Tlc,j ^ kj V T 


X 



(74) 


It is now important to realise that r^(r) is strongly- 
peaked around the condensate edges, hence it is sufficient 
to consider the scaling of the dimensionless condensate 
density n c j = nc,jP 3 in this region. To obtain a quan¬ 
titative estimate, we approximate the condensate by a 
Thomas-Fermi profile and use the density at a healing 
length from the Thomas-Fermi radius as a reference to 
conclude that n c j oc i as k varies. The net result is that 


1/(wt£?) oc l/£ oc 0J c* 0b hence a = 1/2. 

It is now straight forward to see that 1/{ujtc) does not 
scale with (. because of the product n c jn Ct k in Eq. (65), 
hence a = 0. For the miscible Rb-K mixture in Fig. |5| 
this prediction breaks down because the assumption that 
T<c(r) is localised in space is not longer valid; see the left 
middle panel of Fig. [3] 

While the scaled hydrodynamic parameter 77/770 ap¬ 
pears to be small for the exchange collisional process C 12 
when compared to other C 12 and C 22 processes, the ac¬ 
tual numerical values of 77 are in fact large, hence C 12 
remains a dominant process in the situation considered 
by Fig. [5} 


E. Trap geometry 

Besides variation of the isotropic trap frequency, a 
highly-relevant possibility, both experimentally and the¬ 
oretically, is the variation of the trap aspect ratio At ra p = 
w z /o;_l, where lo z (w_l) is the axial (radial) trap frequency, 
so as to probe the physics of reduced dimensionality. If 
Atrap < 1 , we have a cigar-shaped condensate that can 


be used to study, e.g. solitons inn nu and solitonic 
vortices [142] ; if instead, A tr ap > 1, the condensate cloud 
is pancake-shaped and it is commonly used to study vor¬ 
tices HMSJiHa. 

In the following, we choose a reference frequency ui = 
2n x 20Hz and fix lo z (lo±) = lo for cigar (pancake) clouds 
and consider miscible mixtures with two different aspect 
ratios for each trap geometry: 1/A tra p = 05,10 for a 
cigar shaped cloud (column 1,2 in Figs. [ 6 ] and [f| and 
Atrap = 05 ,10 for a pancake shaped cloud (column 4,5 in 
Figs. [ 6 ] an d 0 - We have also checked that our conclusions 
are applicable to immiscible mixtures (not shown). 

Figures [ 6 ] shows the variation of the hydrodynamic 
parameters of 8 'Rb- 41 K mixtures for various collisional 
processes as a function of T/T®, where Tj? is chosen to 
be the critical temperature of the non-interacting single¬ 
component Bose gas for the convenience of comparison. 
For each binary mixture, in going from left to right, the 
trap geometry changes from a quasi-ID geometry to an 
isotropic system, then to a quasi-2D geometry. As the 
geometry changes, all hydrodynamic parameters increase 
like A 0 p for the cigar-shaped cloud, and like yj Atrap for 
the pancake-shaped cloud, meaning that the collisional 
time scale is mainly determined by the tighter trap fre¬ 
quency. This can be understood as a consequence that 
atoms are confined to a smaller region in space with a 
tighter trap, hence an increased probability of collisions. 
Despite the change in the collisional time scales, the rela¬ 
tive magnitudes of 77 when compared among the different 
processes remain roughly unchanged. In other words, if 
the 87 Rb intraspecies scattering is the dominant C 22 pro¬ 
cess in an isotropic trap (third column of top panels in 
Fig. it remains so even if we tighten either the ra¬ 
dial or axial trap frequency. It also means that, the C 12 
collisional process (bottom green solid lines) remains the 
dominant interspecies collisional process (others are in¬ 
dicated by black dots and dash-dots) when the aspect 
ratio is changed. Similar observations can be made on 
87 Rb- 8t> Rb mixtures (Fig. 0 - However, a comparison 
between Fig. [ 6 ] and Fig. [T] reveals another important fea¬ 
ture: the relative magnitudes of 77 , when compared be¬ 
tween intraspecies and interspecies collisional processes, 
are largely determined by the scattering lengths, as we 
have already noted in our analysis of the spatial colli¬ 
sional rates (see the end of Sec. |VB ). For this reason, 
85 Rb intraspecies collisions (blue dashes in the right pan¬ 
els of Fig. [7| dominates both the Ci 2 and C 2 2 processes. 

We would like to caution the reader that our numeri¬ 
cal results are obtained by assuming the non-condensed 
atoms behave like classical particles in three-dimensional 
space. This is certainly not true when the confining 
trap frequency is much larger than the thermal en¬ 
ergy, Juo _l ksT (cigar-shaped cloud) or hco z fcsT 
(pancaked-shaped cloud). Our results here serve only as 
a general guide in changing trap geometry and should not 
be extended to extremely large or small aspect ratios. 
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FIG. 6. (Color online) Variation of the hydrodynamic parameters 1 /(cot) with increasing aspect ratio Atrap = w 3 /ui (left to 
right) for miscible 87 Rb- 41 K mixtures. The scattering lengths are the same as Fig. pi The reference frequency is a) = 2-7T x 20Hz 
and the axial (radial) trap frequency oj z (wi) is equal to cj for quasi-ID (quasi-2D) system. For the ease of comparison, T° is 
chosen to be the critical temperature for the non-interacting single-component trapped gas. 
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FIG. 7. (Color online) Same as Fig. [6] but computed for 87 Rb- 85 Rb mixtures. The scattering lengths are the same as Fig. [ 4 ] 


F. Temperature-dependence of Fc 


The final question that we would like to address with 
our equilibrium simulations concerns the temperature de¬ 
pendence of the exchange collisional process C 12 . At 
equilibrium, the spatially-resolved collision rate (65) can 
be simplified to 


Tc = cr/cj 


M 


kj 


Ulkj 


n c ,k n c ,jv ' r /'(/' + 1), (75) 


where v' r = \j2([Ul — — [U£ - p%])/Mkj is the effec¬ 

tive relative speed and the phase-space distribution is 


/' = 


ex P[(2F77 + ~ vi)/(k B T)} - 1 


(j 


a,b) (76) 


with p chosen such that Eq. ( [76| holds for both j = a, b. 

Since the C 12 collisional process is localised around 
the condensate edge, where U — /i 3 c tends to be small, 
it is tempting to assume that the spatially-resolved col- 
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87 Rb- 41 K 


87 Rb- 85 Rb 



FIG. 8. (Color online) The estimation parameter £ = (p 2 /2rrij + U r ] — pi) /(fcflT) (top) and the exchange collision rate Tc ( |75| ) 
evaluated with the exact phase-space distribution ( |76| ) (bottom black solid) or a phase-space distribution expanded to leading 
order in 1/T (bottom red dashed) as a function of distance r from the trap centre. Shaded regions at the top panels give the 
regime where the high-temperature expansion is inapplicable (£ > 1). Simulation parameters of s 'Rb- 41 K (left) and 87 Rb- 85 Rb 
(right) mixtures are the same as Figs. [3] and [4] at a temperature T = 21nK~ 0.5T C . 


lision rate Ft(r) lies within the high-temperature re¬ 
gion (p 2 /2mj + U/ — p J c -C k B T). In this case, Tay¬ 
lor expansion of /' in orders of 1 /(k B T) then leads to a 
temperature-dependent rate 

r c ~ A(ciT + c 2 T 2 ), (77) 

where A, ci and c 2 are factors to be determined (see 
Appendix |E[). 

In Fig. [STwe show T® (bottom solid) and the estima¬ 
tion parameter (top) 

C = (p 2 / 2m, + Ui - pi)/(k B T) (78) 

as a function of distance r from the trap centre. We 
have 87 Rb- 41 K (left) and 87 Rb- 85 Rb (right) mixtures in 
an isotropic harmonic trap (frequency w = 2-7T x 20Hz) 
at a temperature T = 21nKa; 0.5T C . For better compar¬ 
ison and a somewhat indirect link to approaches based 
on “classical” distribution functions [WHTTU1 . we also 
expand the phase-space distribution to leading order in 
1/T, i.e. /' « k B T/(^~ + U[ — p J c ) and plot the ap¬ 
proximated Tc as dashed lines in the bottom panels. For 
each mixture, the left and right columns show data for the 
miscible and immiscible phases respectively. The bottom 
panels clearly indicate that the high-temperature expan¬ 
sion is valid for the immiscible but not for the miscible 
phase. In the latter case, this is mainly because T® ex¬ 
tends over a broader region in space. 


VI. TWO-FLUID HYDRODYNAMICS 


In this section we derive the hydrodynamic equations 
for the normal components of the multi-component sys¬ 
tem. One of the key theoretical successes of the single 
component ZNG theory is its agreement with the hy¬ 
drodynamic (Landau-Khalatnikov) equations represent¬ 
ing the interaction of the condensed and non-condensed 
components of the system [124 . 145j . This coupling of the 
two fluids has recently been explored for a two component 
Bose system mu, as well as for spin-orbit coupled ther¬ 
mal Bose gases m ■ The hydrodynamic equations for 
the condensate are obtained using the Madelung trans¬ 
formation along with Eq. & yielding 

+ V • (n C j v C j) = - (rf 2 + rj* + r^ 7 ), (79) 


where v c j = (defines the superfluid velocity of 
component j. Equation (79) above defines the continuity 
equation. The Euler equation for component j takes the 
form 


dt 


Vv„- = -VmZ, 


(80) 


where p 3 c = H 2 (V 2 y/n c j) / (2my/n Ct j) + U/ is the non¬ 
equilibrium chemical potential for component j. To ob¬ 
tain the corresponding equations for the non-condensate, 
we take moments with respect to the powers of the mo¬ 
mentum p, i.e. f dp p n (where n = 0, 1, 2) of the kinetic 
equation, Eq. ( |l0| . This leads to a set of three coupled 
nonlinear equations for the non-condensate that can be 
used to describe the limit where collisions dominate, i.e. 
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the hydrodynamic regime. The first of these describes 
the conservation of mass of component j, 


W) = 



dp P(p,r,t), 


(90) 


d /h + V • (n,v nj ) = Tf 2 + + T%, (81) 


respectively. Using eqs. (79) and (81) one finds in turn 
that 


where the velocity of component j is defined as v n j (r, t ), 
while the non-condensate density is fij{ r, t). The velocity 
of component j of the normal fluid is 


f dp p / J (r,p,f) 
J (2tt h) 3 rrij hj(r,t ) 


(82) 


and hj(r,t) has been defined previously by Eq. ©■ 
The corresponding conservation law for the momentum 
of component j (or Navier-Stokes equation) appears as 


d 


+ 


dt 

dp 


d 


m j™j ( is + v n,i ' V ) v n#i j = - n. 


m 

' dx,. 


(27 iky 


P - 771 jV n 


c(i + c% + cfn, (83) 


and (y, p) subscripts refer to spatial components here and 
in what follows. Meanwhile, the conservation law for the 
energy density of component j, ej (r, t) is written 


dej 


+ 


V • (ejV nJ ) = -V • Qj - D 
■l 2 


dt ' ' ^ nj 

dp (p — rrijV, 
(27 tH) 3 2mj 


p 

I *",] 1 fiv ,3 


°12 


C 


kj 


c 


kj 


(84) 


§ i (N c , j (t) + N j (t))= 0. 


(91) 


The conservation of energy and momentum is slightly dif¬ 
ferent, as these quantities are conserved over both com¬ 
ponents, which is reflected in the fact that the delta func¬ 
tions appearing throughout the kinetic theory depend 
generally on both the j and k indices. 

Equations <[79|-((80|) for the condensates and (l8l|), (l83|) 


and (841 for the non-condensates are a direct generaliza¬ 


tion of the equivalent expressions (see Ref. [1451 ) for the 
single-component case, albeit now for a dynamically cou¬ 
pled system comprising two condensates and two ther¬ 
mal clouds. Due to their inherent complexity, the study 
of these coupled equations lies beyond the scope of the 
present work. 

It is anticipated in future works that the hydrodynamic 
equations will yield novel physics, particularly for the 
case of the full Landau-Khalatnikov theory, where the 
entropy of the normal component will cause additional ef¬ 
fects not present in single component thermal Bose gases. 


VII. COMPARISON OF SCHEMES 


The symmetric rate-of-strain tensor appearing in Eq. 
(84) above is defined as 


( r > t) — 2 ^ 


<9v. 


n 


dx u 


+ 


dv 


nv,J 


dx., 


(85) 


and obeys J Dwj_= V • v n .j. The set of equations, 
Eqs. (81), ([83]) and ( |84| introduce the three important lo¬ 
cal hydrodynamic quantities; the stress tensor P Mt/J (r, t), 
the heat current Q ;/ (r, t) and the energy density e, (r, t) 
defined respectively for component j as [40] 


P v v , j— m j 


f dp 

\ Pm 1 

' Pu 

) (2irh) 3 

v nu, 7 

V ni/, 7 

L 771 3 J 


mj f dp 

P 

2 

p 

' 2 J (2 7Th) 3 [ 

v n ,1 

m.j J 


V n,7 

."U J 

f dp 1 



2 

J (2tt h) 3 2mj 

p - mjV nj 

P- 


f j , (86) 
f j , (87) 
( 88 ) 


Note that there is no dependence on the thermal-thermal 
collisional rates T 2 2 in these equations, a consequence of 
the conservation of number, energy and momentum of 
the thermal atoms. To demonstrate this, we can calculate 
the (time-dependent) number of condensate and thermal 
atoms in component j as 



(89) 


Here, we present a brief overview of the different ap¬ 
proaches used to describe the dynamical evolution of cou¬ 
pled multi-component condensates. In our scheme, we 
have, as usual, separated the slowly evolving degrees of 
freedom from those evolving on more rapid time scales. 
In particular, having identified the condensate and (di¬ 
agonal) non-condensate density as the slowly varying 
“mean-field” quantities, we obtained a coupled kinetic 
theory for both condensate and non-condensate that in¬ 
cludes all relevant scattering channels. An important 
point here, similar to some extent to Ref. m is that 
we have treated diagonal elements of the normal pair 
density (corresponds to thermal population) separately 
from their off-diagonal elements (corresponds to “coher¬ 
ences” ). Such a rationale is valid only in the absence of 
optical couplings, where collisions are expected to be the 
dominant process, which is certainly the case for mix¬ 
tures of different species, such as 87 Rb- 85 Rb and 87 Rb- 
41 K considered here, where no interconversion is possible. 
Indeed, for the physically-distinct case of condensates 
with internal spin degrees of freedom, such an adiabatic 
treatment has to be modified in order to self-consistently 
account also for the (internal) coherent coupling of spin 
degrees of freedom 11141 [TT51 . 

There are several other complementary approaches 
to tackling the coupled dynamics of multi-component 
condensates, including the number-conserving ap¬ 
proach m, the stochastic Gross-Pitaevskii formal- 
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ism mmnMi nm . as well as classical field m and 
truncated Wigner |104H106j treatments. 

The starting point of the number-conserving method 
is the Penrose-Onsager criterion, in which the single¬ 
particle density matrix is written in terms of quantum 
field operators. One then expands the field operators 
with the Beliaev decomposition explicitly, maintaining 
their operator form and condensate-noncondensate or¬ 
thogonality, in order to identify the small parameter of 
the theory, namely the (number) ratio of non-condensate 
to condensate population. This in turn allows a set of 
dynamical equations to be extracted which describe the 
condensate through a set of generalized Gross-Pitaevskii 
equations (GGP) and the non-condensate with modified 
Bogoliubov-de Gennes equations. However, the purely- 
dynamical single-component case 1149 1 has yet to be ex¬ 
tended to the multi-component setting. 

In contrast, stochastic treatments of condensate 
dynamics are appropriate for describing the high- 
temperature regime, where the number of atoms in the 
non-condensate fraction is large compared to those in 
the condensate, such that the energetic parameter £ de¬ 
fined by (781 is less than one. Interestingly, the nature 
of multi-component systems means that novel transport 
processes, for example spin-changing collisions will con¬ 
tribute to both damping processes as well as noise for 
these systems. Both the stochastic (projected) GPE 
and the closely-related classical field m and trun¬ 
cated Wigner 111) llilOlij approaches have an explicit cut¬ 
off in energy, with all higher-lying (pure thermal) atoms 
treated stationary. In contrast to these, the ZNG for¬ 
malism explicitly treats all modes dynamically and self- 
consistently, but it does not include the effect of fluctua¬ 
tions of the phase of the non-condensate atoms, which in 
turn limits its applicability to temperatures not too close 
to the transition. 

It is interesting to note that both the multi-component 
number conserving work of m as well as the stochastic 
treatment of nnzi have independently identified a con¬ 
densate to non-condensate “exchange” collisional event, 
physically analogous to that presented in this work by 
C* 2 . In the former case such terms appear as in the 
present work as off-diagonal normal pair averages of fluc¬ 
tuation operators in the corresponding generalized Gross- 
Piteavskii equation, with such averages however defined 
in terms of number-conserving operators, for example see 
Eq. (66a) in [111] . In the stochastic treatment, the ex¬ 
change process enters as a novel scattering amplitude rep¬ 
resenting an extension to the “scattering term” of single¬ 
component stochastic GPE, which however involves mut- 
limode classical field populations, as opposed to those of 
the single-mode condensates arising within our current 
treatment, see Eqs. (73) and (74) in [ 107] . Clearly, each 
of these non-equilibrium theories is quite different in ori¬ 
gin, assumption and applicability, yet all three nonethe¬ 
less demonstrate universal aspects of quantum transport 
theory in low temperature multi-component Bose gases. 

One should also explicitly comment on the link of 


our present approach to the earlier works of Nikuni et 
al. [114 . fl5j . Both works were aimed at discussing spinor 
condensates, for which the Hamiltonian contains addi¬ 
tional terms explicitly maintaining the coupling between 
the two (spin-|) or three (spin-1) different states of the 
system. Clearly, in the case of explicit coupling between 
different states, which enables inter-conversions (i.e. par¬ 
ticles from one state transferred to another state through 
external coupling), our fundamental assumption of treat¬ 
ing the off-diagonal normal pair averages (S^Sk) (j ^ k) 
differently from (S^Sj) could break down. This in turn 
would imply that we should revisit our “slowly-varying 
master” variables, and include {S^Sk) at the same footing 
as (S}jdj) and \4>j\ 2 . This is precisely what has been shown 
in Refs. (1141 IllS] , and would presumably apply to any 
single-species multi-component condensates \F, mp), in 
identical F and different mj? states, under the presence 
of internal or external coupling between those states. Due 
to the added complexity of dealing with an off-diagonal 
Wigner distribution operators, such a model has however 
never been numerically analysed, remaining nonetheless 
an impressive analytical work for such immensely com¬ 
plicated systems. 

In contrast, our present multi-component treatment is 
intended for mixtures of two systems of different species, 
such as a the case of 8 'Rb- 85 Rb and 8 'Rb- 41 K we have 
analysed here, with the full dynamical treatment pend¬ 
ing. 


VIII. CONCLUSIONS 

We have demonstrated how a finite temperature the¬ 
ory describing the out-of-equilibrium dynamics of binary 
Bose gases can be derived using quantum kinetic the¬ 
ory. In particular, it was demonstrated how dissipative 
Gross-Pitaevskii equations for the binary system feature 
three types of collision exchange terms with the non¬ 
condensate atoms of the multi-component system. The 
non-condensates on the other hand are modeled with 
quantum Boltzmann equations coupled to collision inte¬ 
grals which describe atomic scattering between the con¬ 
densate and non-condensate atoms. It was shown in de¬ 
tail how all of these transport processes are derived, in¬ 
cluding the important “exchange” term as well as the 
triplet correlation functions, which are explicitly com¬ 
puted for the multi-component system. 

We also presented results from numerical simulations 
of various binary condensate systems in both miscible 
and immiscible regimes. In particular the condensate 
fractions were estimated for mixtures consisting of 8 ' Rb- 
41 K and 87 Rb- 85 Rb, which showed a slight deviation from 
the estimations based on the single-component Bose gas 
due to mean-field effects. The role of time scales on colli¬ 
sions was elaborated on, in particular the collisionless and 
hydrodynamic regimes were studied through the hydro- 
dynamic parameters of the various collisional processes. 
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Our numerical results demonstrate the interesting 
possibility to access different hydrodynamic regimes. 
For example, thermal-thermal collisions and thermal- 
condensate collisions can occur on comparable or vastly 
distinct time scales through scaling the trap frequency 
and the temperature by the same factor, because the 
various collisional integrals obey different scaling laws. 
On the other hand, the intraspecies collisions can domi¬ 
nate over the interspecies collisions because of the large 
intraspecies scattering lengths, as demonstrated by the 
miscible 87 Rb- 85 Rb mixture, such that Feshbach reso¬ 
nances within and between components could prove use¬ 
ful in investigating different regimes. It is also possible 
to increase the hydrodynamic parameters of all processes 
by changing the trap geometry. However, it is important 
to note that this only changes overall values, rather than 
the relative estimations of different collisional processes, 
which remain unaffected with the C 12 process remain¬ 
ing as the dominant one. We have also investigated the 
extent to which a commonly-implemented phase-space 
expansion to leading order in 1/T is appropriate, finding 
it to be a rather poor estimation in the miscible case. 

The possibility of controlling the hydrodynamicity al¬ 


lows us to explore the interplay of the various collisional 
processes. The hydrodynamic equations developed in 
Sec. [VT] shows the added complexity due to the extra 
scattering channels present in the binary system. With 
the presence of eight collisional processes (three C 22 pro¬ 
cesses, four C 12 processes and a C 12 process) in a bi¬ 
nary system compared to two collisional processes (a C 22 
process and a C 12 process) in a single-component Bose 
gas, it is not a priori clear how an out-of-equilibrium bi¬ 
nary mixture would relax to its final state and at what 
time scales, especially if several collisional processes were 
close to the hydrodynamic regime and were competing 
with each other. A definitive answer requires careful nu¬ 
merical simulations of the full non-equilibrium dynamics, 
which is a subject under our active investigation. 
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Appendix A: Perturbing Hamiltonian 


The purpose of this appendix is to detail the steps required to arrive at the perturbing Hamiltonian. After applying 
Wick’s theorem to H 3 and H 4 , it can be shown that the non-quadratic terms appearing in the full system Hamiltonian 
H can be recast in the form H 3 —>• 5H\ and H 4 —> SH 2 A SHq respectively, where 
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where hkj = (Sj^k) and fhkj = {djSk}- Here Eqs. (Ala)-(Ale) are generated by inserting Eqs. (16a) and (16b) into 
Eqs. (15d) and ( |15e| ). The term given by Eq. (Ala) constitutes a mean-field shift to the chemical potential of the 
system. However, as we follow the usual prescription of keeping interaction effects within the chemical potential to 
linear order in the scattering length | |124| . these terms need not be considered any further. With the definitions of 
Eq. (Ala)-(Ale I, one can show that the resulting perturbing Hamiltonian 


H\t) = H - Hmf 


(A2) 


can be broken into contributions comprising different numbers of fluctuation operators. Using Eq. the perturbing 
Hamiltonian can be expressed in the form [2 Di 1124 ] 


H’{t) = H[(t) A H' 2 (t) A H' 3 (t) A 
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where 


H[{t) = -5H, = - dr\ 22 
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The full expression for the Fourier transformed perturbing Hamiltonian becomes 

H\t) = H[(t) + H' 2 {t) + H' 3 (t) + H' 4 (t), 
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which are used to derive the collision integrals in the body of the text. Next, each term in Eq. (A5) is decomposed 
into an intra- H' n j{t) and inter- H' nk -{t) component contribution so that for n= 1—4 one has the decomposition 
H' n {t) = H' n j(t) + H’ nk -{t). Then using the Fourier expansion of the fluctuation operator 5j(r,to) one obtains 
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The expressions given above by Eqs. (A6a)-(A6h) allow us to derive the collisional integrals in the body of the text. 
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Appendix B: Anomalous Pair Correlation Functions 

In this appendix the anomalous pair correlation functions of the form (SjSj) and (SkSj) appearing in Eq. |6]) are 
calculated. Following the steps detailed in the body of the paper, one can shown that the pair anomalous terms take 
the form 


(& k Sj ) — I TT <jl : j 0 k 0j ( ^ 




S { £J C+ £ C- £ 1-4) S P 3 
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By taking the continuum limit, one finds that Eq. (Bl) becomes the collisional integral 
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valid for both j=k and j ^ k This integral describes a collision process whereby a pair of condensate particles collide 
to give two thermal particles and it’s inverse process. Such terms have not been included in our theory as they violate 
energy conservation. For this reason, these terms are conventionally dropped due to the Popov approximation m- 
Figure [9] shows schematically the different kinetic scattering processes represented by Eq. B2 Squares represent 
condensate and circles represent non-condensate (thermal) atoms, while the blue and red colors represent the a and 
b components respectively. 


Pi P2 P2 Pi 

( i ) m' s + 

P3 P P P3 


(hSj) 


FIG. 9. (Color online) The two diagrams represent kinetic: (i) intra -component ( S a 5 a ); (h) inter-component: (5 a Sb) pair 
anomalous scattering terms of Eq. (B21. Corresponding diagrams for the b component are obtained by interchanging the 
(open) red and (closed) blue colors appearing above. 


Appendix C: Evaluation of Non-Equilibrium averages 

Here, the form of various non-equilibrium averages used in the intermediate steps of deriving the full collisional 
theory are listed. We begin by writing out in full the triplet anomalous averages which are used in the definitions of 
the condensate growth terms R^ and R k ' 3 as 


<Wj>(D=-i / rft , (5 t (t / ,to)[5 t (t,t , )44^5(t,t , ) J ^ )fcj (*)]^(i , ,io)), 
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and 
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*0 

with similar expressions for j = k. Meanwhile, the condensate exchange terms IR.^- 7 are calculated as (for j ^ k) 

t 

$&>=-J /(C3) 


This expression, based on symmetry-breaking did not appear in previous works for spinor gases HUGE], but does 
formally appear in the related c-field approach of Ref. m■ 
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Next, we give the form of the quantum Boltzmann contributions C\ 2 and Cyj required to give the final form of the 
associated collision integrals. These are 


C{ J 2 = —Tr / 5(t,to)[/' 7 (r,P^o),-ff3,iW] I 


and 
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Using the Fourier transform of the Wigner operator 
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along with H 3 (Eq. (A6f|), it can be shown that the right hand side of Eq. (C5| can be written as 
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The expression given above contains two distinct types of three-field correlation functions, both of which can be 
computed using the definition of the non-equilibrium average along with the relevant contribution from the perturbing 
Hamiltonian, Eq. (A6f). One can in particular show that 


'J,P2 < U.P3^,P4) — + )^( £ c+ e p2 £ p 3 £ p 4 )^Pf + P2, 


P3+P4 


/l(/I+i)(/i+i)-(/Ei)/l/i 


(C8) 


and the corresponding expression for (aj p4 aj P3 «j, P2 ) can be obtained by taking the Hermitian conjugate of Eq. (C8). 
To simplify the first summation over the center of mass momentum q appearing in Eq. (C7), we can expand the sum 
and note that terms with q ^ 0 make zero overall contribution. After which, upon inserting Eq. (C8) into Eq. (C7) 


results in the final expression given by Eq. (45) of Sec. IV C 1 

Meanwhile, the exchange collisional term cf J 2 is computed as (for j k) 
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Then, two quadratic correlation functions are required in order to write down a final expression for the collision 
integral, Eq. (CIO). As before we can use the expression for the non-equilibrium average, Eq. (251 along with the 
relevant part of the perturbing Hamiltonian, Eq. (A6d). This gives 
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Inserting Eq. (Cll) into (CIO) above allows us to write down the discrete form of the exchange collision integral, Eq. 
(49) of Sec. IV C 2 of the text. 

The final non-equilibrium average that is required are those terms describing scattering exclusively between thermal 
atoms, C ° 2 2 and C 2 %, which is given by the expressions 


C 22 = ^Trp(M 
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Proceeding, we can write down an expression for C^.j using (Eq. (A6h|) which takes the form 
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Again to simplify Eq. (C14), we require the quartic correlation function comprised of equal numbers of creation and 
annihilation operators. This is computed using the definition given by Eq. (251 along with H ’ 4 k -(t) as defined by Eq. 
(A6h). This gives 
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Interchanging the dummy summation variables over momentum in Eq. (C14), and inserting the expression given by 
Eq. (C15) leads to Eq. (531 of Sec. IV C 3 of the main text. 

To obtain all of the collision integrals appearing in the body of the text, one also requires the following approxima¬ 
tions 
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Finally, in evaluating the integrals over time appearing in all of the collisional integrals, the following identity has 
been used 
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where V (...) represents the Cauchy principle value in Eq. (C19) above, and is conventially dropped [111 ;. 115 . 


Appendix D: Transformation between lab frame and centre-of-mass frame 


In the numerical evaluation of the collisional rates (60), (61) and (62), it is useful to transform the momenta from 
the lab frame to the centre-of-mass frame. To this end, we define the center-of-mass and relative momenta, (P, p r ) 
and (P', p(,), such that 
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The two Dirac delta functions then enforce P = P', |p r | = jp(,| because of energy and momentum conservation. 
The collision rates between non-condensed atoms (for both k = j and k ^ j) is 
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where akj = (1 + 5kj)^a\j is the cross section, Vi and V2 are the initial velocities of atoms j and k respectively, 12 
specifies the solid angle of the final relative velocity V 4 — V 3 . 

For collisions between condensate and non-condensate atoms, we first consider the process (for both k = j and 
k j), the “out” collision rate that represents the scattering of a non-condensed atom from a condensate to produce 
two non-condensed atoms is given by 
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We again use (D1) to obtain 
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where n° ut = y | v c j — v 2 | 2 — 2([/„ — n 3 c )/mkj is the relative speed of the initial states corrected to take into account 
the local conservation of energy. 

The reverse process, where two non-condensed atoms collide such that one of them goes into a condensate, is given 
by the “in” rate as 
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In order to reduce Eq. (D6) into a useful form for Monte Carlo sampling as well as dynamical simulations, we follow 
the approach of Jackson and Zaremba |127| and arrive at 
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where v™ = V4 — v J c is the velocity of thermal atom j relative to the local condensate velocity while the second integral 
is a two-dimensional integral over the velocity vector v which is in a plane normal to vj, n . The velocity of the other 
incoming thermal atom v 3 is then given by 
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and the outgoing velocity of the thermal atom is given by 
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Note that we follow and drop the cubic term 72/3/4 in numerical simulations because it cancels exactly between 
the “in” and “out” rates. 

For the exchange collisions (65), instead of the usual center-of-mass transformation (Dl), we use an alternative 
transformation 
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The exchange collision rate is therefore 
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where M k 1 = m k 1 — m • 1 plays the role of an effective reduced mass while the effective relative speed is 
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Appendix E: Equilibrium evaluation of 


In general the collision integrals cannot be evaluated analytically. However, the exchange rate, Eq. (Dll) can be 

2 _ _ _ 

calculated in the limit Po(^~ +Uj! — <C 1. Using the momentum transformations given by (DlOa) and (DlOb), 

the scattering rate FJ ut becomes 
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Finally, the ‘relative’ momentum is defined as p ont = \Jp 2 + 2A4 k j(gjjn c j — g kk n c ,k)- Equation (El) shows the 


temperature dependence of r£; 7 8 9 10 11 12 , while the constants ci and c 2 defined by Eq. (E2) and (E31 give respectively the 
general form of the coefficients of T 2 and T. 
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